#'
#' Find all of the individals in the annulus of given focal individuals
#'
#' @param data a population or community object
#' @param fx,fy the location of focal individuals
#' @param type the shape of quadart, circle or rectangle
#' @param rRange the minimun and maximun of cricle ring quadrat
#' @param qsize the length and width of rectangle quadrat
#' @param info what kinds of information should be return, "index": index of individual only, "dist": distance only, "both"
#' @param keepraw whether save the raw fx, fy and data within the result.
#' 
#' @details
#' The fixed-radius nearest neighbor algorithm was inplemented here.
#' 
#' 
#' @return
#' A frnn object contains the indexes of neighbors and distances to neighbors
#' 
#' @author
#' Guochun Shen, Changjiang Gou
#' 
#' @examples
#' data(BCI)
#' focusx=BCI$x
#' focusy=BCI$y
#' focusxy=matrix(c(focusx,focusy),nrow=total_abundance(BCI),ncol=2)
#' 
#' system.time(re1<-frnn(BCI,focusx,focusy,rRange=c(0,5)))
#' system.time(re1<-frnn(BCI,focusx,focusy,rRange=c(0,5),info="index"))
#' library(spdep)
#' system.time(re2<-dnearneigh(focusxy,d1=0,d2=5))
#' 
#' checkpoint=sample(1:total_abundance(BCI),1)
#' id1=sort(re1[[checkpoint]])
#' id2=sort(re2[[checkpoint]])
#' 
#' id1
#' id2
#' 
#' plot(x=BCI$x[id2],y=BCI$y[id2])
#' points(x=BCI$x[id1],y=BCI$y[id1],pch=19,col=2)
#' 
#' plotdim=attr(BCI,"plotdim")
#' fx=runif(100,0,plotdim[1])
#' fy=runif(100,0,plotdim[2])
#' 
#' re3=frnn(BCI,fx,fy,rRange=c(0,5))
#' 


#' @export
frnn=function(data,fx,fy,rRange=c(0,10),qsize=c(0,0),type="circle",info="both",keepraw=FALSE){
  plotdim=attr(data,"plotdim")
  
  if(any(fx<0 | fy<0 | fx>plotdim[1] | fy>plotdim[2]))
    stop("some focus points located outside of the plot")
  
  #due to the problem of code in c++, we need to enlarge the plotdim a little bit
  plotdim=plotdim+min(plotdim)*0.01
  
  xycoli=match(c("x","y"),colnames(data))
  xy=as.matrix(data[,xycoli])
  if(type=="circle"){
    typei=1
  }else{
    typei=2
  }
  if(info=="both"){
    infoi=0
  }else if(info=="index"){
    infoi=1
  }else{
    infoi=2
  }
  
  result=inner_sample(xy,plotdim[1],plotdim[2],fx,fy,rRange[1],rRange[2],qsize[1],qsize[2],typei,infoi)
  
  class(result)=c("frnn",class(result))
  if(keepraw){
    attr(result,"focus_xy")=data.frame(x=fx,y=fy)
    attr(result,"neighbors")=data
  }
  return(result)
}

#' @export
plot.frnn=function(frnn,wait=FALSE){
  
  focus_xy=attr(frnn,"focus_xy")
  neighbors=attr(frnn,"neighbors")
  
  
  if(is.null(focus_xy) | is.null(neighbors)){
    stop("No raw spatial distribution data were stored in the given object, please set 'keepraw=TRUE' in frnn")
  }
  
  if(!(c("index") %in% names(frnn[[1]])))
    stop("Index information is needed from frnn, please set info='both' or info='index' in frnn")
  require(ggplot2)
  point_data=data.frame(x=c(focus_xy[,1],neighbors$x), y=c(focus_xy[,2],neighbors$y),
                      type=rep(c("Focus","Neighbors"),times=c(nrow(focus_xy),nrow(neighbors))))
  
  
  nindex=unlist(lapply(frnn,function(x) x$index))
  arrow_data=matrix(nrow=length(nindex),ncol=4)
  arrow_data[,3]=neighbors$x[nindex]
  arrow_data[,4]=neighbors$y[nindex]
  nn=lapply(frnn,function(x) length(x$index))
  arrow_data[,1]=rep(focus_xy[,1],nn)
  arrow_data[,2]=rep(focus_xy[,2],nn)
  arrow_data=as.data.frame(arrow_data)
  colnames(arrow_data)=c("x0","y0","x1","y1")
  
  if(nrow(arrow_data)>1000){
    warning("Too much  (>1000) arrows to show, might be very slow!")
    if(!wait) stop()
  }
  
  fig=ggplot()+
    geom_point(aes(x=x,y=y,col=type),data=point_data)+
    geom_segment(aes(x=x0,y=y0,xend=x1,yend=y1),arrow = arrow(length=unit(0.03,"npc")),data=arrow_data)+
    theme_bw()
  return(fig)
}


# the R version
frnnR=function(data,fx,fy,rRange,qsize,info="both",type="circle"){
 plotdim=attr(data,"plotdim")
 N=total_abundance(data)
 Nindex=1:N
 Nfocal=length(fx)
 if(type=="circle"){
   r1=rRange[1]
   r2=rRange[2]
 }else{
   r1=qsize[1]
   r2=qsize[2]
 }
 
 grid_size=floor(r2*2)
 xbreaks=unique(c(seq(0,plotdim[1],by=grid_size),plotdim[1]))
 ybreaks=unique(c(seq(0,plotdim[2],by=grid_size),plotdim[2]))
 nxi=length(xbreaks)-1
 nyi=length(ybreaks)-1
 ngrid=nxi*nyi
 
 #put focus points into the grid
 fxi=as.numeric(cut(fx,breaks=xbreaks,include.lowest=TRUE))
 fyi=as.numeric(cut(fy,breaks=xbreaks,include.lowest=TRUE))
 fgridi=fxi+(fyi-1)*nxi
 
 
 #put neighbor points into the grid
 xi=as.numeric(cut(data$x,breaks=xbreaks,include.lowest=TRUE))
 yi=as.numeric(cut(data$y,breaks=xbreaks,include.lowest=TRUE))
 gridi=xi+(yi-1)*nxi
 
 #summarize neighbor points into different grid
 points_grid=lapply(1:ngrid,function(x) Nindex[gridi==x] )
 
 result=list()
 for(i in 1:Nfocal){
   #calculate neigbor grids
   posnxi=c(fxi[i],fxi[i]+1,fxi[i]-1)
   posnyi=c(fyi[i],fyi[i]+1,fyi[i]-1)
   posnxi=posnxi[posnxi>0 & posnxi<=nxi]
   posnyi=posnyi[posnyi>0 & posnyi<=nyi]
   posgridi=rep(posnxi,times=length(posnyi))+(rep(posnyi,each=length(posnxi))-1)*nxi
   
   #get the indexes of individuals in the neighbor grid and focal grid
   posindex=unlist(points_grid[posgridi])
   #calculate spatial distsances from focal individual to its neighbors
   ndist=sqrt((data$x[posindex]-fx[i])^2+(data$y[posindex]-fy[i])^2)
   
   #select individuals satisfy conditions
   if(type=="circle"){
     
     sel= ndist>r1 & ndist<=r2 
   }else{
     sel=(data$x[posindex]-fx[i])<=r1 &  (data$y[posindex]-fy[i])<=r2
   }
   
   if(info=="both"){
     selindex=posindex[sel]
     seldist=ndist[sel]
     result[[i]]=list(index=selindex,dist=seldist) 
   }else if(info=="index"){
     selindex=posindex[sel]
     result[[i]]=selindex
   }else{
     seldist=ndist[sel]
     result[[i]]=seldist
   }
   
 }
 class(result)=c("frnn",class(result))
 return(result)
}

#' @useDynLib scp
#' @importFrom Rcpp sourceCpp
NULL